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FOREWORD 


The  Investigation  reported  herein  is  part  of  the  research  project  at 
The  Ohio  State  University,  Columbus,  Ohio  supported  by  the  Air  Force  Of¬ 
fice  of  Scientific  Research  Grant  83-00-55.  Lt.  Col.  John  J.  Allen  was 
the  Program  Manager  at  the  commencement  of  the  project.  Lt.  Col.  Law¬ 
rence  D.  Hokanson  was  the  Program  Manager  from  July  1,  1983.  The  re¬ 
search  project  was  started  Feburary  1,  1983  and  is  continuing.  The 
present  report  documents  part  of  the  work  done  up  to  January  31,  1984. 
At  The  Ohio  State  University,  the  project  is  supervised  by  Dr.  Ranbir 
S.  Sandhu,  Professor,  Department  of  Civil  Engineering.  The  computer  pro¬ 
gram  modification  and  analyses  reported  herein  were  started  by  Dr.  Baher 
L.  Aboustit,  Post-doctoral  Reseach  Associate  and  completed  by  Mr,  Soon- 
jo  Hong,  Graduate  student,  who  also  prepared  the  documentation  on  the 
Ohio  State  University  Computer  System.  The  Instruction  and  Research 
Computer  Center  at  The  Ohio  State  University  provided  the  computational 
facilities. 


ABSTRACT 


Numerical  performance  of  Ghaboussi  and  Wilson's  isoparametric  bili¬ 
near  quadrilateral  element,  for  analysis  of  quasi-static  flow  of  an  in¬ 
compressible  fluid  through  a  linear  elastic  saturated  porous  soil,  is 
compared  with  that  of  Sandhu's  composite  element  in  which  the  displace¬ 
ment  has  biquadratic  interpolation.  Application  of  both  procedures  to 
solution  of  Terzaghi's  one- dimensional  consolidation  and  Gibson's  prob¬ 
lem  of  plane  strain  consolidation  of  the  half-space  under  a  strip  load 
shows  that  Ghaboussi  and  Wilson's  procedure  gives  results  almost  identi¬ 
cal  to  those  from  the  higher  order  element  but  is  significantly  more  ec¬ 
onomical  to  use. 
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SECTION  I 


INTRODUCTION 

Since  the  first  application  of  the  finite  element  method  to  analysis 
of  seepage  in  elastic  media  (1,2)*  considerable  progress  has  been 
made  in  the  theoretical  formulation  as  well  as  computational  procedures. 
Recent  advances  include  variational  formulations  admitting  limited 
smoothness  of  finite  element  bases  (3,4),  experimentation  with  several 
different  spatial  interpolation  schemes  and  investigation  of  various 
temporal  approximation  methods  (5-15).  The  finite  element  method  has 
been  applied  to  soils  exhibiting  secondary  compression  (8,11,13),  nonli¬ 
near  soil  behavior  (13,14-18),  and  to  finite  deformation  (16,18).  Devel¬ 
opments  in  solution  procedures  include  use  of  Laplace  transforms 
(10,11,19),  automatic  selection  of  the  time-step  size  (13)  and  the  use 
of  single  variable  formulations  (20,21). 

In  spatial  discretization,  Sandhu  (1)  proposed  that  the  order  of 
terms  appearing  in  a  convolution  product  in  the  variational  principle  be 
the  same.  This  produced  the  "composite"  element  in  which  the  order  of 
polynomial  interpolation  for  displacements  was  higher  than  that  for 
fluid  pressures.  The  composite  element  first  proposed  by  Sandhu  (1,2), 
and  used  by  Hwang  (22)  and  others,  was  the  "63"  element  with  quadratic 
interpolation  for  displacements  and  linear  interpolation  for  fluid  pres- 
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sures  over  triangular  regions.  Later,  Sandhu  (7,23)  introduced  the  "84" 
element  which  had  eight  point  biquadratic  interpolation  for  displace¬ 
ments  and  was  a  four  point  isoparametric  quadrilateral  for  fluid  pres¬ 
sures.  This  element  was  also  used  by  Runesson  (13).  Buchmaier  (24)  ex¬ 
perimented  with  five,  six  and  seven  point  quadrilaterals  for  fluid 
pressure  as  transition  elements  near  loaded  surfaces. 

Several  spatial  interpolation  schemes,  besides  the  composite  ele¬ 
ments,  have  been  tried  by  various  investigators.  Ghaboussi  (5)  used 
four  point  isoparametric  quadrilaterals  for  both  the  fields.  However,  an 
additional  incompatible  mode  was  included  in  the  displacement  approxima¬ 
tion.  This  element  has  the  economy  while  the  additional  "local"  mode 
gives  it  the  character  of  a  "higher  order"  scheme.  Smith's  (25)  formula¬ 
tion  was  similar  to  Ghaboussi 's  except  that  no  incompatible  modes  were 
used.  Prevost  (9)  proposed  cautious  use  of  "reduced  integration"  in  con¬ 
duction  with  Smith's  "44"  element.  Yokoo  (26),  Booker  (19),  and  Vermeer 
(27)  used  triangular  elements  with  linear  interpolation  for  both  the 
displacement  and  the  fluid  pressure  fields.  Other  interpolation  schemes, 
based  on  use  of  quadrilaterals  built  up  from  linear  four  or  five  point 
triangles  were  tried  by  Sandhu  (8). 

All  investigators  have  generally  reported  success  with  whatever 
scheme  they  used.  Comparative  studies  of  different  elements  are  rare. 
Some  comparisons  of  nemerical  performances  were  attemped  by  Sandhu 
(7,8).  In  evaluating  various  candidate  schemes,  Sandhu  proposed  that  an 
acceptable  method  meet  the  following  requirements  in  addition  to  effi¬ 
ciency  and  accuracy. 

i.  The  interpolation  scheme  must  conform  with  the  assumptions 
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regarding  continuity  and  differentiability  used  in  setting  up  the 
governing  variational  formulation. 

ii.  It  should  be  possible  to  generate  the  "undrained"  solution,  i.e. 

the  state  of  fluid  pressures  and  displacement  at  time  t  =  0+. 
iii.  For  sufficiently  small  time  steps,  the  scheme  should  be  insensi¬ 
tive  to  the  choice  of  the  time-step  size. 

Elements  "63"  and  "84"  satisfy  these  requirements.  However,  the  compos¬ 
ite  elements  are  too  expensive  to  be  used  in  large  problems.  This  has 
discouraged  extension  of  the  analysis  to  three-dimensions,  and  to  nonli¬ 
near  and  dynamic  problems.  The  4-4  element  is  more  economical  but  was 
found  by  Sandhu  (8)  to  have  oscillatory  errors.  The  "63"  element  has 
been  widely  used  because  it  was  the  first  one  to  be  introduced  and  gave 
satisfactory  results  in  most  cases.  However,  the  8-4  element  gives  re¬ 
sults  almost  identical  to  those  from  the  "63"  element  but  is  more  eco¬ 
nomical  as  it  requires  fewer  nodal  points  and  has  less  band-width.  Gha- 
boussi's  element  (5)  was  not  included  in  this  comparison. 

The  purpose  of  the  present  investigation  was  to  implement  Ghaboussi's 
formulation  (we  refer  to  it  as  the  6-4  element  because  of  the  additional 
internal  modes  in  the  displacement  representation)  and  to  compare  its 
performance  with  that  of  the  8-4  element.  Both  the  elements  were  used  to 
solve  Terzaghi's  problem  of  one-dimensional  consolidation  and  Gibson's 
problem  of  a  half  space  under  a  strip  load.  To  simulate  consolidation  of 
the  half  space  by  a  finite  domain  model,  an  approximation  to  the  condi¬ 
tions  at  infinity  is  required.  Three  alternative  assumptions  regarding 
the  "cut  off"  boundary  were  considered. 


SECTION  II  summarizes  the  governing  field  equations  following 
Ghaboussi  (5).  The  finite  element  development  Is  given  In  SECTION  III. 
SECTION  IV  describes  the  example  problems.  Results  of  the  analysis  are 
given  In  SECTION  V. 


SECTION  II 


EQUATIONS  GOVERNING  LINEAR  ELASTIC  SOIL  CONSOLIDATION 


Assuming  pore  water  to  be  compressible,  the  equations  of  force  equi¬ 
librium  of  elementary  volumes  and  mass  continuity,  over  the  spatial  re¬ 
gion  of  interest  R,  may  be  written  in  standard  indicial  notation  as, 
[Appendix  A], 


(1) 


U,j(7r,  *  *  “j.j  ' 


(2) 


where  u. ,  f.,  E.,..,  K..,  denote  the  cartesian  components,  respectively, 

1  1  ixllj  1J 

of  the  displacement  vector,  the  body  force  vector  per  unit  mass,  the  iso¬ 
thermal  elasticity  tensor  and  the  permeability  tensor. pis  the  mass  den¬ 
sity  of  the  saturated  soil  and  p^  that  of  water,  n  is  the  pore  water 
pressure,  a  is  the  solid  compressibility  and  M  is  a  measure  of  fluid 
compressibility.  With  these  field  equations  we  associate  the  following 
boundary  conditions; 
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Here,  t^,  are  components  of  the  traction  and  fluid  flux  vectors  asso¬ 
ciated  with  surfaces  embedded  in  the  closure  of  R.  x  ■  are  components  of 

^  J 

the  total  stress  tensor.  $2^  are  complementary  subsets  of  the  bound¬ 
ary  of  the  spatial  region  of  interest  and  so  are  S^,  S^.  Even  though  the 
equations  given  above  apply  to  compressible  fluids,  the  applications  re¬ 
ported  herein  assumed  incompressible  fluid  i.e.  M  oo  .  The  Initial 
conditions  for  the  problem  are 


0)  =  u.(0)  on  R 
7r(x  0)  =7r(0) 

J 


on  R 


SECTION  III 


FINITE  ELEMENT  FORMULATION 

Discretization  of  the  governing  function  for  the  two-field  formula¬ 
tion  followed  by  application  of  the  variational  principle  (5)  leads  to 
the  following  matrix  equation. 


Fk  k  1 

uu  pu 

U(tj) 

■  b  o' 

Pi' 

< 

’+ 

1 

<  [-OAtK  -C  ] 
Lpu  pp  pp:i 

ir(ti)_ 

_  -Kp„  C-d-ajAtKpptCpp] 

.P2. 

where  (^Q.t^)  is  the  single  time  step  of  interest,  and; 

At  »  tj-tg 

f|u(to)]=  vectors  of  nodal  point  values  of  the  components  of 
the  displacement  at  time  respectively 

|7r(ti)|  ,|7r(to)j  =  vectors  of  nodal  point  values  of  the  pore  water  pre¬ 
ssure  at  time  respectively. 

^Pj}=  the  vector  of  nodal  point  loads  including  applied  nodal  loads, 
boundary  tractions,  body  forces,  initial  stresses  and  effect  of 
displacement  constraints. 

{^2)-  vector  of  nodal  point  fluxes  including  applied  nodal  fluxes, 
boundary  fluxes,  body  force  effects  and  effects  of  specified 
pore  water  pressures. 

[•^uu^“  the  spatial  "stiffness  matrix"  for  the  elastic  soil. 

CKpp]=  the  spatial  "flow  matrix"  for  the  compressible  fluid  and  At=  1. 
a  =  the  coefficient  characterizing  single-step  temporal  discretl- 
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zation. 

[•^pu]=  the  coupling  matrix  representing  the  influence  of  pore  pressure 
in  the  force  equilibrium  equation. 

the  coupling  matrix  representing  the  influence  of  soil  volume 
change  upon  the  nodal  point  fluid  pressure. 
tCpp]=  the  spatial  fluid  compressibility  matrix. 

The  matrix  depends  upon  the  interpolation  scheme  for  displacements 

3nd  CCpp],  [Kpp]  depend  upon  the  interpolation  scheme  for  the  pore-water 
pressures.  The  coupling  matrix  [Kp^]  involves  spatial  interpolation  for 
both  the  field  variables.  The  temporal  discretization  for  the  single 
step  scheme  is  reflected  in  the  value  of  the  coefficient  a.  For  linear 
interpolation  a  =  0.5. 

Equation  (7)  includes  the  "natural"  boundary  conditions  expressed  by 
Equation  (4)  and  (6).  Equation  (3)  and  (5)  are  satisfied  by  explicitly 
requiring  u^  =  u.jOn  S^^and  77’='^  on  S3.  Development  of  the  vectors  and 
matrices  appearing  in  Equation  (7)  is  given  in  Appendix  B.  As  the  fluid 
is  assumed  to  be  incompressible  in  both  the  example  problems,  the  matrix 
Cpp  is  zero.  It  is  included  in  the  theoretical  discussion  for  complete¬ 
ness.  The  computer  program  can  handle  compressible  as  well  as  incompres¬ 
sible  fluids. 


•  »  •  •  n  •  ^ 

•  *  m  ^ 


r  t  •  . 

V-VJ 


SECTION  IV 


EXAMPLE  PROBLEMS 

An  existing  computer  program,  initially  developed  by  Dr.  J.K.  Lee  of 
the  Department  of  Engineering  Mechanics.  The  Ohio  State  University,  was 
modified  to  include  plane  strain  consolidation  using  both  the  8-4  and 
the  6-4  elements.  The  code  was  used  to  solve  Terzaghi's  problem  of  one¬ 
dimensional  consolidation  and  Gibson's  problem  of  two-dimensional  con¬ 
solidation.  For  these  problems,  the  theoretical  solutions  are  known  and, 
therefore,  precise  comparison  was  possible. 

For  Terzaghi's  problem,  the  dimensions  of  the  consolidating  soil  col¬ 
umn  and  soil  properties  were  the  same  as  in  Sandhu's  example  (23);  i.e. 
soil  depth  h  =  7,  modulus  of  elasticity  E  =  6000,  Poisson's  ratio  = 
0.4,  coefficient  of  permeability  K  =  4x10  .  Soil  and  water  particles 
were  assumed  to  be  incompressible,  i.e.  o  =  1,  M->oBin  Equation  (7). 
Figure  1  illustrates  the  problem  and  Figure  2  illustrates  the  mesh  used 
in  the  analysis.  The  mesh  consisted  of  9  elements  with  20  nodes  when 
using  the  6-4  element  and  48  nodes  when  using  the  8-4  element.  In  tem¬ 
poral  discretization,  the  coefficient  o  was  given  the  value  of  0.5.  The 
following  temporal  partition  was  considered. 

10  steps  of  At  =  0.000001  over  [0,0.00001] 

10  steps  of  At  =  0.01  over  [0.00001,0.1] 

10  steps  of  At  *  0.1  over  [0.1,  1.1] 

10  steps  of  At  =  10  over  [1.1,  101.1] 

10  steps  of  At  =  100  over  [101.1,  901.1] 


where  unit  is  second 


The  second  problem  was  that  of  a  saturated  layer  of  finite  thickness 
subjected  to  a  strip  loading.  An  analytical  solution  for  this  problem 
was  developed  by  Gibson  et  al.  (28).  The  material  properties  were  as¬ 
sumed  the  same  as  for  the  one- dimensional  problem  except  that  Poisson's 
ratio  was  set  equal  to  zero  in  the  problem.  The  geometry  and  the  finite 
element  mesh  are  shown  in  Figure  3.  The  mesh  consisted  of  88  elements 
with  108  nodes  when  using  the  6-4  element  and  303  nodes  when  using  the 
8-4  element.  Since  the  infinite  domain  had  to  be  modeled  by  a  finite 
one  in  the  numerical  schemes,  the  following  alternative  boundary  condi¬ 
tions  were  specified  at  the  far  boundary  where  the  half-space  was  cut 
off. 

a.  Pore  fluid  pressure  and  horizontal  displacement  prescribed. 

b.  Fluid  flux  and  horizontal  displacement  prescribed. 

c.  Pore  pressure  and  traction  prescribed. 

In  temporal  discretization,  the  factor  o  was  set  equal  to  0.5,  i.e. 
piecewise  linear  variation  in  u^  andTTwas  assumed.  The  following  par¬ 
tition  was  used; 

10  steps  of  At  =  0.1  over  [0,1] 

10  steps  of  At  =  1  over  [1,11] 

10  steps  of  At  =  10  over  [11,111] 

10  steps  of  At  =  100  over  [111,1111] 

10  steps  of  At  =  1000  over  [1111,11111] 


Here,  unit  is  second. 


SECTION  V 


RESULTS  OF  ANALYSIS 

A.  Terzaghi's  Problem  of  One-Dimensional  Cosol idation 

In  analyzing  Terzaghi's  problem.  Figure  1,  CPU  time  on  the  AMDAHL 
470/V8  Computer  was  13.92  seconds  for  8-4  element  while  for  the  6-4  ele¬ 
ment  it  was  10.03.  Table  1  shows  time  settlement  history  for  the  two 
types  of  elements  as  well  as  the  analytical  solution.  The  response  of 
the  two  elements  practically  coincided  throughout  the  time  domain  and 
was  in  good  agreement  with  the  analytical  solution,  except  at  early  time 
i.e.  t  <  0.1.  The  pore-pressure  distribution  generated  by  the  two  ele¬ 
ments  for  different  time  steps  is  shown  in  Figure  4.  The  vertical  axis 
shows  the  ratio  of  depth  below  surface  to  the  total  height  h  of  the  soil 
column  considered.  The  horizontal  axis  gives  the  pore  fluid  pressure  as 
a  fraction  of  the  intensity  of  the  constant  surface  load.  The  plots  are 
for  values  of  time  variable  equal  to  0.000001,  11.1,  201.1,  and  901.1 

corresponding  to  non-dimensional  time  factor  T  =  Ct/h  of  0.1x10  , 

0.01165,  0.2110,  0.9195  where  C  =  2GK  and  K  is  the  coefficient  of  perme¬ 
ability.  The  plots  of  pressure  distribution  history  generated  by  the  two 
schemes  coincide.  At  early  stages,  the  error  in  the  pore  pressure  at 
points  near  the  loaded  surface  is  quite  large  for  both  the  schemes.  This 
is  a  feature  of  the  spatial  interpolation  (8)  used. 


TABLE  1 


Surface  Settlement  History 


Time 

(sec) 

8-4  Element 

6-4  Element 

Exact(Ref.23) 

0.2 

0.51917x10-5 

0.51583x10-5 

0.28147x10-5 

0.6 

0.15881x10-5 

0.15849x10-5 

0.15417x10-5 

1.1 

0.21221x10-4 

0.21189x10-4 

0.20874x10-4 

21.1 

0.91394x10-4 

0.91365x10-4 

0.91423x10-4 

41.1 

0.12815x10-3 

0.12813x10-3 

0.12760x10-3 

81.1 

0.18022x10-3 

0.18019x10-3 

0.17924x10-3 

301.1 

0.34446x10-3 

0.34450x10-3 

0.34205x10-3 

901.1 

0.50352x10-3 

0.50351x10-3 

0.50166x10-3 

Q  TIME=  0. 1000-05 
0  TIME=  0.1110+02 
A  TIME=  0.2010+03 
A  TIME=  0.9010+03 


Figure  4:  Fluid  Pressure  Distribution  at  Various  Time  Stages 
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B.  Gibson's  Problem  of  Plane  Strain  Consolidation  under  a  Strip  Load 

In  analyzing  the  second  problem,  the  CPU  time  on  the  AMDAHL  470/V8 
Computer  for  8-4  element  was  315  seconds  while  for  the  6-4  element  it 
was  33  seconds.  Table  2  shows  vertical  displacement  and  pore  pressure  at 
the  distance  0.084a  below  the  surface  directly  under  the  center  of  the 
loaded  area.  Here,  a  is  half-width  of  the  loaded  strip.  A  non-dimension¬ 
al  measure  of  the  settlement  is  introduced  as  Gu/aq  where  G  is  the  shear 
modulus,  q  the  intensity  of  load  per  unit  width  of  the  strip.  The  fluid 
pressure  is  represented  by  the  fraction xr/q.  The  quantities  are  tabulat- 
ed  for  dimensionless  time  given  by  T=  Ct/hr  Good  agreement  between  the 
results  for  the  two  finite  element  procedures  is  seen  at  all  time  stag¬ 
es.  Table  3  and  Figure  5  show  comparison  between  the  settlement  and  pore 
pressure  history  for  cases  a  and  b,  using  the  6-4  element  under  diffe¬ 
rent  conditions  at  the  "cut  off"  boundary.  From  Table  3  and  Figure  5,  it 
appears  that  the  location  of  the  "far"  boundary  was  selected  far  enough 
so  that,  for  the  case  where  horizontal  displacements  are  assumed  to  van¬ 
ish  at  this  distace,  prescribing  vanishing  fluid  pressure  or  fluid  vel¬ 
ocity  at  this  boundary  had  little  effect  on  the  displacements  and  fluid 
pressures  near  the  center  of  the  loaded  area. 

Comparing  case  c  with  case  a,  i.e.  considering  the  effect  of  traction 
or  displacements  at  the  "far"  boundary  for  prescribed  fluid  pressure,  it 
was  found.  Figure  6,  that  the  settlement  in  the  early  stage  of  loading 
was  significantly  affected.  The  solution  for  case  c  was  in  excellent 
agreement  with  Gibson's  (28)  theoretical  solution. 


TABLE  2 


Vertical  Displacement  and  Pore  Pressure  at  0.084a  below  the  Center  of 
the  Loaded  Area.  Case  a:  Fluid  Pressure  Prescribed  at  the  "far" 

Boundary 


9.6x10’^ 


3.84x10 


0.01056 


0.03936 


0.10656 


0.39456 


0.68256 


1.06656 


Gu/aq 

8-4 

6-4 

0.176136 

0.177012 

0.180666 

0.181698 

0.190212 

0.191256 

0.218784 

0.22014 

0.262584 

0.264534 

0.365376 

0.36753 

0.413892 

0.416008 

0.443682 

0.445914 

8-4 


0.88874 


0.48602 


0.32230 


0.17829 


0.11498 


0.50406 


0.25815 


0.01117 


6- 


0.9002 


0.48610 


0.32164 


0.18030 


0.11526 


0.50271 


0.25778 


0.01117 


3.94656 


0.466842 


0.469134 


0.1174x10 


0.115x10 


I 


TABLE  3 

Comparison  of  Case  a  and  b:  Vertical  Displacement  and  Pore  Pressure  at 
0.084a  below  the  Center  of  the  Loaded  Area 


Ct 

3.6x10-5 

Gu/aq 

TT/q  1 

Case  a 

Case  b 

Case  a 

Case  b 

0.17855 

0.177012 

0.89619 

0.90026 

9.6x10-^ 

0.17926 

0.17769 

0.73221 

0.73562 

0.01056 

0.192738 

0.191256 

0.32014 

0.32164 

0.04896 

0.229296 

0.22779 

0.16341 

0.16434 

0.10656 

0.265962 

0.264534 

0.11456 

0.11526 

0.87456 

0.433548 

0.434202 

0.017178 

0.016863 

1.06656 

0.445332 

0.445914 

0.011439 

0.011163 

SECTION  VI 


CONCLUSIONS 

Ghaboussi  and  Wilson's  6-4  element  and  Sandhu's  8-4  element  were  ap¬ 
plied  to  Terzaghi's  and  Gibson's  problem  for  which  exact  solutions  are 
available.  Results  of  these  limited  tests  show; 

i.  The  6-4  element  gave  a  solution  identical  to  that  given  by  the  8-4 
element  but  with  significant  savings  in  computational  time. 

ii.  The  6-4  element  as  well  as  the  8-4  element  predicted  the  un¬ 
drained  solution  (t  =  0+). 

iii.  At  early  stage  of  loading,  both  elements  gave  unsatisfactory  re¬ 
sults.  Apparently,  special  singularity  elements  (29)  are  required  near 
loaded  drained  surfaces, 

iv.  The  conditions  at  the  "cut  off"  boundary  for  the  finite  domain 
considered  in  the  finite  element  model  must  be  carefully  defined.  The 
solution  near  the  loaded  region  was  not  sensitive  to  whether  the  fluid 
flux  or  fluid  presure  was  prescribed.  However,  prescribed  displacement 
or  tractions  had  significant  influence.  This  is  in  line  with  Sandhu  and 
Schiffman's  experience  (30). 

V.  The  6-4  element  is  distinctly  superior  to  the  4-4  element  in  that 
it  gives  the  solution  at  time  =  0+  and  does  not  have  the  oscillatory  er¬ 
ror  of  the  4-4  element.  At  the  same  time,  it  has  the  economy  of  the  sim¬ 
pler  element. 
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vi.  Using  the  8-4  element  to  generate  the  solution  at  time  t  =  0+, 
numbering  of  the  nodal  points  has  to  be  carefully  ordered  to  avoid  ze¬ 
roes  on  the  diagonal  of  the  matrix.  However,  Ghaboussi  and  Wilson's  6-4 
element  is  free  from  this  defect.  This  is  because,  in  eliminating  the 
additional  degrees  of  freedom,  the  static  condensation  would,  in  gener¬ 
al,  result  in  non-zero  diagonal  quantities. 

vii.  Finally,  it  appears  that  Ghaboussi  and  Wilson's  6-4  element  is  a 
good  candidate  for  consideration  towards  application  in  nonlinear  and 
dynamic  problems  as  well  as  for  extension  to  three-dimensional  applica- 
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APPENDIX  A 


EQUATIONS  OF  CONSOLIDATION 

In  this  appendix,  we  list  the  field  equations  governing  the  quasis¬ 
tatic  deformation  of  soil-water  mixtures,  following  Biot's  theory  (31). 


A.l  PRELIMINARIES 

Let  R  be  the  open  connected  region  occupied  by  the  fluid-solid  mix¬ 
ture.  Let  OR  be  the  boundary  of  R  and  IT  its  closure.  The  domain  of  defi¬ 
nitions  of  all  functions  of  interest  is  the  cartesian  product  R  x  [0,oo)^ 
Here,  [0,oo)  denotes  the  positive  interval  of  time. 


Let  e.  (x,t),  Ui(x,t),  w,-(x,t),  e,'j(x,t),  rTj(£.t),  f^(x,t), 


<l^(x,t)  and  0j(x»t),  in  that  order,  denote  cartesian  components  of  the 
displacement  vector,  the  fluid  displacement  vector,  the  relative  dis¬ 
placement  of  the  fluid  with  respect  to  the  solid,  the  infinitesimal 
strain  tensor,  the  symmetric  Cauchy  effective  stress  tensor  in  the  soil, 
the  total  stress  tensor,  the  body  force  vector,  the  fluid  flux  vector 
and  the  quantity,  conjugate  to  the  fluid  flux  vector,  denoted  by  vector 
0.  Also,  let  7r(x,t)  denote  the  isotropic  fluid  stress  above  a  reference 
state  (throughout,  we  shall  consider  the  atmospheric  pressure  to  be  the 
reference  state).  Let  p  denote  the  bulk  density  of  the  mixture, 
p equals  the  sum  of  the  equivalent  bulk  density  Pj  of  the  soil  and  p^  of 
the  fluid.  Following  Biot  (31),  we  assume  the  stresses  d(x,t)  and 
TTtx.t)  to  act,  respectively,  over  the  solid  and  the  fluid  areas  of  any 
surface  element. 


A. 2  THE  FIELD  EQUATIONS 
A. 2.1  Kinematics 

For  small  strains,  infinitesimal  strain  tensor  is  expressed  in  terms 
of  displacement  components  as 


•ij  =  ‘^2  <“t.j  *  “j.t)  =  “(i.j) 


(A-1) 


Displacement  of  the  fluid  relative  to  the  solid  but  measured  in  terms  of 
volume  per  unit  area  of  the  bulk  medium  is  defined  by 


”i  ^  f(Ui  -  Ui) 


(A-2) 


where  f  is  the  porosity.  Accordingly,  for  uniform  porosity,  the  fluid 
relative  volumetric  strain  is  given  by 


^<^i,i  -  ^i,i) 


(A-3) 


A. 2. 2  Constitutive  Equations 

The  constitutive  equations  have  to  be  written  for  the  stresses  as 
well  as  for  diffusive  resistance. 


A. 2. 2.1  Total  and  Pore  Pressure 

Biot  (31),  assuming  the  existence  of  a  strain  energy  function  for  the 
mixture,  postulated  the  following  constitutive  relations  for  compressi¬ 
ble  fluid  following  through  linearly  elastic  solid. 


Tij  =  Enkieki  +  OTrSij  on  TT  x  [0,oo) 


(A-4) 


M  (a  e  +  ^)  on  R  X  CO.oo)  (A-5) 


where  a,  M  are  material  coefficients. 


A. 2. 2. 2  Diffusive  Resistance 

Diffusive  resistance  in  a  binary  mixture  of  soil  and  water,  for  no 
inertia  effects,  is  defined  (1)  by  the  components 


Vf  s  -JT,  -  P  f, 


(A-6) 


(A-7) 


For  the  linear  theory,  assuming  irrotational  isothermal  flow,  the  diffu¬ 
sive  resistance  is  linearly  proportional  to  the  relative  velocity  of  the 
fluid-solid  constituents.  This  leads  to 


^,i  ■*■^2^1  "  ®i  =  -  =  ^ij”j 


(A-8) 


where  w^  are  components  of  the  nominal  relative  velocity  and  are 
are  components  of  the  flow-resistivity  tensor.  Symmetry  requirements  im- 


(A-9) 


If  the  resistivity  tensor  has  an  inverse  denoted  by  components  K.., 

^  J 

Equation  {A-8)  gives  a  generalization  of  Darcy  law 


Wi  =  qi  = 


S' 


(A-10) 


The  permeability  tensor  is  symmetric  so  that 


(A-11) 


A. 3  BALANCE  LAMS 
A. 3.1  Force  Equilibrium 

Balance  of  forces  on  infinitesimal  volumes,  assuming  small  deforma¬ 
tion,  neglecting  inertia  and  in  the  absence  of  body  couples  is  expressed 
by  equations 


T-  •  •  +  Pf  •  =0 
HJ.J  1 


on  R  X  C0,ae) 


■Jl 


(A-12) 

(A-13) 


A, 3.2  Mass  Continuity 

For  no  chemical  reaction  between  constituents  of  a  binary  mixture, 
the  mass  of  each  within  a  fixed  volume  must  be  conserved.  Referring  to 
material  frame  associated  with  the  soil  skeleton,  for  compressible 
fluid,  continuity  impies  continued  saturation  of  the  soil.  Formally,  we 
write 


K. .( 


TT  .  +p„fj 


2' j  M 


Combining  Equations  (A-5)  and  (A-14)  to  eliminate 


(A-14) 


■O'd  *  l-WV.  K,j(7rj  ‘Pjfj),, 


(A-15) 


i.e.,  for  compressible  fluid  and  solid,  the  rate  of  the  outflow  from  a 
material  volume  of  soil  skeleton  equals  the  rate  of  decrease  of  the  vo¬ 
lume.  Equation  {A-12)  upon  substitution  of  Equation  (A-4)  and  (A-2) 
yields  Equation  (1)  of  SECTION  II.  Equation  (2)  in  that  section  is  a  re¬ 
arrangement  of  Equation  (A-15). 


APPENDIX  B 


GHABOUSSI  AND  WILSON’S  ELEMENT 

In  this  appendix,  we  outline  the  development  of  matrices  for  the  6-4 
element. 

B.l  GEOMETRY 

The  cartesian  coordinates  of  any  arbitrary  point  in  a  quadrilateral 
may  be  written  as 


4 


(B-1) 

1*1 

4 

y  ■  Z»i"i 

i-1 

(B-2) 

where  N^.  are  the  interpolating  functions  and  ,  y^  are  the  coordinates 
of  the  four  nodal  points  of  the  quadrilateral. 

B.2  DISPLACEMENT 

For  an  isoparametric  element,  the  interpolating  functions  for  the 
displacement  and  geometry  are  the  same.  However,  in  this  element  two  ex¬ 
tra  modes  are  used  in  defining  displacements,  i.e. 


3( 


where  u  and  v  are  the  components  of  displacement  vector  along  x  and  y 

direction,  are  the  nodal  point  displacements  or  the  generalized 

coordinates  associated  with  the  additional  modes. 

B.3  PRESSURE 

The  fluid  pressure  at  an  arbitrary  point  of  the  isoparametric  element 
is 


4 

*  Z’^iNi  =  Nff?  (B-5) 

1-1 

where  y  is  the  vector  of  nodal  pressure.  In  Equations  (B-1)  through 
(B-5),  the  ^  i=l,6  are 

Nj  =  1/4  (l-s)(l-t) 

N2  =  1/4  (l+s)(l-t) 

N3  =  1/4  (l+s)(l+t) 

(B-6) 

N4  =  1/4  (l-s)(l+t) 

N5  =  (l-s^) 

Nj  •  (1-t^) 

where  s,  t  are  the  local  coordinates  of  the  element  such  that  the  oppo¬ 
site  pairs  of  sides  of  the  element  are  represented  by  s  =  ±1  and  t  =  +1. 
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B.4  STRAIN-DISPLACEMENT  RELATIONSHIP 


D/Dx  0|  fu 


fiyy  '  =  0  way  ‘ 
2e^y  [o/Oy  a/Dx]  [v 


6 


8/ay  a/ax 


1-1 


6 

,x''i'*'Ni  ,yUi ). 


=  Bp  U 


(B-7) 


where  U  is  the  vector  of  nodal  and  generalized  displacement  and  B  is  an 
appropriate  transformation.  The  volumetric  strain  is  given  by 


^  ®xx  ®yy  ®  ^^i.x'^i  2*^i,y''i  “  -A  ~ 


B.5  PRESSURE  GRADIENT 

Components  of  the  pressure  gradient  are; 


(B-8) 


’®x  h;x 


i]N^.x»'i 


(B-9) 


B.6  EVALUATION  OF  DERIVATIVES  .  AND  N, 

■"  ^  tj 

Equation  (B-6)  gives  the  interpolating  functions  in  terms  of  the  lo¬ 
cal  coordinates  s,t.  To  evaluate  derivatives  with  respect  to  global 
coordinates  x,y,  we  use  the  relationship 


=  1/J 


by/5t  -Sy/bs 
-bx/6t  9x/0s 


(B-10) 


4  4. 

where  .s^i  )(§'i  .t^i )  -  (gi  .s^i  Xgi  .t’^i  )•  Equation  (B-10)  is 

obtained  by  inverting  the  expressions  of  N.  and  N.  .  in  terms  of  N.  ^ 

and  N.  by  the  chain  rule  of  differentiation. 

'  »y 


B.7  ELEMENT  MATRICES 

The  element  matrices  contributing  to  the  system  described  by  Equation 
(7)  are  defined  as  follows. 


M  -I 

'pp 

•I  -I 

"up  =// V'  •'  ■  "pu' 

'pp  -/'/'is-  ^ 


(B-11) 

(B-12) 

{B-13) 

(B-14) 


Here  D,  k  are  matrices  describing  the  elastic  properties  and  the  perme 
ability,  respectively,  of  the  porous  material. 


B.8  LOAD  VECTORS 


The  vector  on  the  right  hand  side  of  Equation  (7)  are 

given  by 

+  [Vi^\  (B-15) 

IP2}  =  aAt  M2(ti)  +  (1-a)  AtM2(to)  -  a  At  M4(ti) 


-(IHI)  At  M4(t^j)  (B-16) 

where 

fC  Ny  [Pf]  J  dsdt  (B-17) 

{”2}  = f  \  k  \Pf)  J  dsdt  (B-18) 

^”3)  =  /*  Nu  Nt  (B-15) 

•'s 

(M  =  jT  N„  10}  dS4  (B-20) 


where  N^,  are  defined  in  Equations  (B-3),  (B-4)  and  {B-5),  and  f  is 
assumed  to  be  constant  over  each  element.  N^,  N^,  N^,  J^represent  inter¬ 
polating  functions  for  boundary  tractions,  boundary  flux,  displacements 
of  the  solid  and  fluid  pressures,  respectively.  The  interpolating  func¬ 
tions  iSj^,  for  the  consol idating  region  and  for  its  boundary  will  be 
different. 


B.9  ELIMINATION  OF  LOCAL  MODES 

For  the  6-4  element,  the  discretized  equations  can  be  written  for  a 
time  step  (^Q.tj)  in  the  form. 


[K  ] 

uu-* 

u{tj) 

■  i 

«u 

[K 

pu-* 

-[0At[K|,p]-[Cpp3 

irCtj) 

•  4 

R,r 

L  y 

where  {R^}  =  iMi)  +  (M3)  (B-22) 

+  (l-a)At[M2(tQ)}  -aAt(M4(ti))  -  (l-a)At{M4(to)}  (B-23) 


To  eliminate  the  four  degrees  of  freedom  associated  with  the  two  incom¬ 
patible  modes  for  each  displacement  component,  the  consolidation  Equa¬ 
tion  (B-21)  can  be  partitioned  as; 

*-'^uull^  •^'^uul2^ 

^'^uu21^  I^'^uu22^ 

/"pul^^  tSu2J^  -(fiAt[Kpp].[Cpp]) 

where  the  internal  degrees  of  freedom  are  denoted  by  u^.  Solving  Equat¬ 
ion  (B-24)2  for  U2, 

“2  *  >^0022  M  *^02  ■  '^uu2l“l  •  '^du2’^] 


(B-25) 


Substituting  for  u^  In  Equation  (B-24)j  and  collecting  terms. 


*-'^uull  ■  *^0012  '^uu22’^ 

*  *-'^pul  ■  *^0012  *^0022  ^  '^pu2^  {?r(tj)} 

“  ^"^ul  ■  *^0012  *^0022  ^  %2} 

Similarly,  substituting  for  in  Equation  (8-24)3  and  collecting  terms 

*^0022  ^  '^uu21^ 

■  "  Sp  "  W'  W2''  Su2^  Wh)} 

=  W'  W2‘^«u2i 

Equations  (B-26)  and  (B-27)  can  be  rewritten  as 


C  C]  I 

_  V  Sp  .  ^ 


(B-28 


where 


*-'^uu^  “  “^uuH  "  '^uul2  *^0022  ^  *^0021 


(B-29 


^*^pu^  “  ’^pu  "  *^uul2  *^0022  *^pu2 

f^'^pp]  *  -aAtKpp  -  Kpy2  Kyy22  ^  l^pu 2  ’  Sp 


(B-3C 


(B-31 


*  ^^ul  "  '^uul2  *^0022  *^02 


(B-32 


^R,r3  “  ■  '^py2^  *^0022  ^  ^^02 


(B-33 


It  is  worth  noting  that  even  if  the  fluid  is  incompressible,  i.e.  Cpp  is 
null,  at  At=0,  Kpp  is,  in  general,  nonzero.  This  makes  solution  of  the 
undrained  problem  possible. 
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